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Abstract 

The paper is devoted to the effects of superconducting pairing in small metallic grains. It 
turns out that at strong superconducting coupling and in the limit of large Thouless conductance 
one can explicitly determine the low energy spectrum of the problem. We start with the strong 
coupling limit and develop a systematic expansion in powers of the inverse coupling constant for 
the many-particle spectrum of the system. The strong coupling expansion is based on the formal 
exact solution of the Richardson model and converges for realistic values of the coupling constant. 
We use this expansion to study the low energy excitations of the system, in particular energy and 
spin gaps in the many-body spectrum. 

1 Introduction 

Since mid 1990's, when Ralph, Black, and Tinkham succeeded in resolving the discrete excitation 
spectrum of nanoscale superconducting metallic grains pQ, there has been considerable effort to describe 
theoretically superconducting correlations in such grains (see e.g. j2] for a review). However, very 
few explicit analytical results relevant for the low energy physics of superconducting grains have 
been obtained, since, in contrast to bulk materials, the discreetness of single electron levels plays an 
important role. In this paper we address this problem in the regime of well developed superconducting 
correlations. 

The electron-electron interactions in weakly disordered grains with negligible spin-orbit interaction 
are described by a simple Hamiltonian [3| 

H un [ Y _ = H BC s - JS(S + 1) (1) 

N 

Hbcs = e i c L c i<? ~ Xd J2 c li c h c jl c jl ( 2 ) 

where e; are single electron energy levels, d is the mean level spacing, cJ CT and q ct are creation and 
annihilation operators for an electron on level i, S and N are the total spin and number of levels 
respectively. There are only two sample-dependent coupling constants: A and J that correspond to 
superconducting correlations and spin-exchange interactions respectively. Throughout the present 
paper, for the sake of brevity, we consider only the less trivial case of ferromagnetic exchange, J > 0. 

Although Hamiltonian (|TJ is integrable [3J |5] and solvable by Bethe's Ansatz, the exact solution 
|Hj yields a complicated set of coupled polynomial equations (see Eq. (jSJ) bellow). As a consequence, 
very few explicit results have been derived and most studies resorted to numerics based on the exact 
solution. The purpose of the present paper is to remedy this situation and to build a simple and 
intuitive picture of the low energy physics of isolated grains in the superconducting phase. 
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It is well known that physical observables of a superconductor are nonanalytic in the coupling 
constant A at A = 0. On the other hand, the opposite limit of large A turns out to be regular and 
relatively simple. Here we use the exact solution to obtain an explicit expansion in powers of 1/A for 
the ground state and low lying excitation energies. 

We will distinguish between two types of excitations: ones that preserve the number of Cooper 
pairs (the number of doubly occupied orbitals) and ones that do not. Only the latter excitations are 
capable of carrying nonzero spin. It turns out that for J = to the lowest order in 1/ A both types of 
excitations are gaped with the same gap XNd. We compute explicitly the two gaps to the next nonzero 
order in 1/A and find the gap for pair-breaking excitations to be larger. The difference between the 
two gaps turns out to be of the order of d 2 / A, where d is the mean single particle level spacing and A 
is the BCS energy gap, i.e. the difference vanishes in the thermo dynamical limit. We were not able 
to determine the convergence criteria for the strong coupling expansion exactly, however we present 
evidence that the expansion converges up to realistic values of A between A c i ~ 1 and A C 2 ~ l/V. 

Hamiltonian (0) was studied extensively in 1960's in the context of pair correlations in nuclear 
matter (see e.g. [7j). A straightforward but important observation was that singly occupied levels 
do not participate in pair scattering [Hj. Hence, the labels of these levels are good quantum numbers 
and their contribution to the total energy is only through the kinetic and the spin-exchange terms 
in Due to this "blocking effect" the problem of diagonalizing the full Hamiltonian reduces to 
finding the spectrum of the BCS Hamiltonian (J2J) on the subspace of either empty or doubly occupied - 
"unblocked" orbitals. The latter problem turns out to be solvable jH] by Bethe's Ansatz. The spectrum 
is obtained from the following set of algebraic equations for unknown parameters Ef. 

1 m n n i 

Xd j=1 Ei — Ej k=1 Ei — 2efc 

where m is the number of pairs and n is the number of unblocked orbitals e^. Bethe's Ansatz equations 
(jSJ) for the BCS Hamiltonian (j2J) are commonly referred to as Richardson's equations. The eigenvalues 
of the full Hamiltonian are known to be related to Richardson parameters, E^, via 

m 

E = Y,Ei + Y,£B-JS(S + l) (4) 

i=l B 

where J2b € b is a sum over singly occupied - "blocked" orbitals and S is the total spin of blocked 
orbitals (i.e. the total spin of the system). 

BCS results [9 j for the energy gap, condensation energy, excitation spectrum, etc. are recovered 
from exact solution (jSJ) in the thermodynamical limit ^U] • The proper limit is obtained by taking the 
number of levels, N, to infinity, so that Nd — > 2D = const, m = n/2 = N/2, where D is an ultraviolet 
cutoff usually identified with Debye energy. In particular, for equally spaced levels e i: the energy gap 
A and the ground state energy in the thermodynamical limit are 

A < A ) = ^iWA) *fr CS = -^»lA < 5 > 

Since the BCS Hamiltonian (J2J) contains only three energy scales: D, A, and d, there are only two 
independent dimensionless parameters: N, and A. The perturbation theory in small A breaks down in 



2 



the superconducting state as is already suggested by BCS formulas (J3J). Thus, it is natural to consider 
the opposite limit of large A and treat the kinetic term in Hamiltonian ((!} as a perturbation. 

The paper is organized as follows. In Section 2 we consider the limit A — > oo, which is the zeroth 
order of our expansion. In this limit one can determine the spectrum straightforwardly by representing 
the BCS Hamiltonian (j2J) in terms of Anderson pseudospin operators [llj . In particular, one finds that 
at J = excitations with nonzero spin to the lowest order in 1/A have the same gap (the spin gap) 
as spinless excitations. Next, we rederive the same results from Richardson's equations and also 
show that in the limit A — > oo the roots of Richardson's equations are zeroes of Laguerre polynomials. 

In Section 3 Bethe's Ansatz equations © are used to expand the ground state and low-lying 
excitation energies in series in 1/A. We write down several lowest orders explicitly and give recurrence 
relations that relate the kth order term to preceding terms. These relations can be used to readily 
expand up to any reasonably high order in 1/A. Finally, we compute the spin gap to the next nontrivial 
order in 1/A and demonstrate that at J = the first excited state always have zero spin. 



2 The strong coupling limit 

In this section we analyze the lowest order of the strong coupling expansion. 

As the strength of the coupling constant A increases, the spectrum of the BCS Hamiltonian (J2J) 
undergoes dramatic changes as compared to the spectrum of noninter acting Hamiltonian Hbcs{^ — 0). 
First, there is a region of small A where the superconducting coupling causes only small perturbations 
in the electronic system. This region shrinks to zero in the thermodynamical limit and is roughly 
determined by the condition A(A) < d jTTj, where A(A) is given by (j3J). For larger A the perturbation 
theory in A breaks down ^3] and strong superconducting correlations develop in the system. A 
representative energy level diagram is shown on Fig 1. In the crossover regime the spectrum displays 
numerous level crossings which reflect the break down of perturbation theory in A. The fact that the 
crossings occur for random single electron levels e^, i.e. in the absence of any spatial symmetry, is a 
characteristic feature of quantum integrability [T4] . 

The lowest order of the strong coupling expansion is obtained by neglecting the kinetic energy term 
in the BCS Hamiltonian (j2J). This limit can in principle be realized in a grain of an ideal regular shape 
|15j . In this case the single electron levels are highly degenerate and if the energy distance between 
degenerate many-body levels is much larger than Xd, only the partially filled Fermi level is relevant. 
Then, the kinetic term in (j2J) is simply a constant proportional to the total number of particles and 
can be set to zero. 

An efficient way to obtain the spectrum of Hamiltonian (0) in the strong coupling limit is by repre- 
senting the interaction term in the BCS Hamiltonian in terms of Anderson pseudospin- 1/2 operators 

m 

K , = + Kt = W)t = c , tC(i (6) 

The pseudospin is defined only on unblocked levels, where it has all properties of spin-1/2, i.e. proper 
commutation relations and definite value of Kf = 3/4. 

The interaction term in the BCS Hamiltonian (j2J) takes a simple form in terms of K{. 



H™ cs = -Xd K + K- = -Xd 



K(K + l)-(K z ) 2 + K z (7) 
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Figure 1: Results of exact numerical diagonalization. Energies of BCS Hamiltonian (J2J) for m = 4 pairs 
and n = 8 unblocked single particle levels versus coupling constant A. All energies are measured 
in units of the mean level spacing d. The single particle levels 6j are computer generated random 
numbers. As the strength of the coupling A increases, the levels coalesce into narrow well separated 
rays (bands). The width of these bands vanishes in the limit A — > oo (see Eq. (|47|) and the discussion 
around it). Slopes of the rays and the number of states in each ray are given by Eq. (JTJ [T31 I12J1 . The 
ground state is nondegenerate, while the first group of excited states contains n — 1 = 7 states. Note 
also the level crossings for A ~ 1 (see the insert on the above graph). 
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where K = J2iKi is the total pseudospin of the unblocked levels. The z-projection of the total 
pseudospin according to (jOJ is K z = m — n/2, where m and n are the total number of pairs and 
unblocked (either doubly occupied or empty) levels respectively. It is simple to check that replacing 
a doubly occupied level with two singly occupied ones does not affect the difference m — n/2. As a 
result, 

n N 

K* = m-~ = M-— (8) 
2 2 v ; 

where M is the maximum possible number of pairs and N is the total number of levels respectively. 
Hence, the last two terms in (J7J) yield a constant independent of the number of blocked levels. This 
constant can be set to zero by an overall shift of all energies. Therefore, the full Hamiltonian (0) in 
the strong coupling limit is 

H umv. = ~^K(K + 1) - JS(S + 1) (9) 
Since there are n pseudospin-l/2s, the total pseudospin K takes values between \K Z \ and n/2, 

->K>\m--l (10) 

while the total spin S ranges from (1/2) to M — m (M — m + 1/2) for even (odd) total number of 
electrons. For the sake of brevity, let us from now on consider only the case of even total number of 
electrons. Then, the sum of the total spin and pseudospin is constrained by 

N 

K + S<- (11) 

The degree of degeneracy D(K, S, n) of each level is [TH] 

n(Kq] nKgg+l) (N-n)\(2S+l) 

1 ' } (f +# + 1)1(2-10! (*=z + S + l)\(Z=z-S)\ 1 ] 

The ground state of Hamiltonian © has the maximum possible pseudospin, K = N/2, and minimal 
possible spin, S — 0, provided that Xd > J (recall that we consider only positive values of the exchange 
coupling J). 

There are two ways to create an elementary excitation. First, one can decrease the total pseudospin 
K while keeping the total number of pairs M unchanged. The second type of excitations corresponds 
to breaking pairs and blocking some of the single electron levels. These excitations can contribute to 
the total spin of the grain S. They also affect the pseudospin since its maximal value K max = n/2 
is determined by the number of unblocked levels. The lowest-lying excitations correspond to K — 
N/2 — 1, which can be achieved both with and without breaking a single Cooper pair. Therefore, 
we find from Q that the pair-conserving excitations are separated by a gap A Da j r = NXd while 
pair-breaking excitations can lower their energy by having nonzero spin S. Since the maximum value 
of S for two unpaired electrons is S = 1, we get A g pj n = NXd — 2 J . 

In the opposite case J > Xd, K = and the total spin has the maximum possible value S = M in 
the ground state, i.e. J = Xd is the threshold of Stoner instability in the strong coupling limit. 

The above results can be obtained directly from exact solution (J3J). Moreover, individual parameters 
Ei can also be determined and, since eigenstates of the BCS Hamiltonian (J2J are given in terms of Ei 
(see 0), this can be used to calculate various correlation functions in the strong coupling limit. 
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The value of the total pseudospin K turns out to be related to the number, r, of those roots of 
equations (jSJ) which diverge in the limit A — > oo (see bellow). To the lowest order in 1/A we can neglect 
single electron levels 6j in Eqs. (jUJ) for these roots 

1 + J2 2 - n ' i - 1 r (13) 
\d . =1 Ei — Ej Ei 

where n' = n + 2r — 2m and summation excludes j = i. For the remaining m — r roots we have 

n y 

^2— — = i = r + l,...,m-r (14) 

k=1 J^i — 



Multiplying each equation in ()T3j) by E{ and adding all Eqs. fl3j) . we obtain the eigenenergies of the 
BCS Hamiltonian (J2J) for n unblocked levels and m pairs 

E = -Xd r(n - 2m + r + 1) (15) 

Comparing this to (J7|) and (JHJ), we find the relationship between r and K 

r = K + m-n/2 (16) 

Since the total pseudospin, K, is constrained by (jlOjl. the number, r, of diverging Richardson param- 
eters, is also constrained 

2m — n < r < m ii n < 2m 

(17) 

0<r<m ifn> 2m 

Bellow in this Section we show that Eqs. (JT5j) have a unique solution. As a result, the degeneracy 
of energy levels ()12j) is equal to the number of solutions of Eqs. (JT4*j) for the remaining E^ This number 
can be computed [T71HB] directly from (JT4*j) and indeed coincides with (fT^j) . 

Finally, Eqs. (|13j) can be solved to determine parameters Ei to the lowest order in 1/A (see also 
[TTH I2TI]). To this end it is convenient to introduce a polynomial f(x) of order r with zeroes at 
x = Xi = Ei/(Xd) 

f{x) = f[(x-Xi), (18) 



i=l 



Using 



one can rewrite Eqs. (fTHj) as 

F(xi) = where F(x) = x/'(a;) - x/"(x) + n'f'(x) (19) 

Since F(x) and f(x) are two polynomials of the same degree r with the same roots x«, they are 
proportional to each other. The coefficient of proportionality is the ratio of coefficients at x r and, 
according to ()19|). is equal to r. Therefore, F(x) = rf(x), or equivalently 

xf" - (x + n')f + r/ = (20) 
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The only polynomial solution to this equation is the Laguerre polynomial L r 1 n ' . Thus, to the order 
A the nonvanishing roots of Richardson's equations <jBJ) in the strong coupling limit are determined by 

L r 1 ~ n> (aJ) = n ' = n + 2r - 2m ( 21 ) 

where r is the number of nonvanishing roots to the order A. This number and the total pseudospin are 
related by (|TE|). The ground state has r = m, the first degenerate group of excited states corresponds 
to r — m — 1, etc. The constraint r > 2m — n in (JT7jl follows from the requirement that the roots of 
(|2*Tjl be nonvanishing Moreover, it can be shown using conditions (j!7|) that all Richardson 
parameters E t are complex for even values of r, while for odd r there is a single real (negative) root. 
The fact that the roots of (|13|) are generally complex was also noted in on the basis of numerical 
solution of Richardson's equations. 



3 The Strong Coupling Expansion 

Now we turn to the expansion in powers of 1/A around the strong coupling limit. The evolution of 
energy levels with A can be viewed as a motion of one-dimensional particles whose positions are the 
energies of the BCS Hamiltonian (J2J) (see e.g. [ZHED- Then, single electron levels determine the 
initial conditions at A = 0. As the coupling A increases beyond the crossover between the weakly 
perturbed Fermi gas and the regime of strong superconducting correlations, the particles gradually 
loose the memory of their initial positions and eventually the spectrum becomes independent of tj. 
In this limit, the excited levels coalesce into highly degenerate rays with a universal slope (see Fig. 1 
and Eq. ((15])) . In the strong coupling expansion the system of one-dimensional particles evolves from 
larger to smaller A. One expects this evolution to be nonsingular until we come close to the level 
crossings (see the beginning of the previous section), i.e. the crossover region, where both expansions 
in A and in 1/A break down. 

A quantitative estimate of the convergence of 1/A expansion can be obtained by considering various 
limiting cases. In the thermodynamical limit the ground state energy is given by BCS expression (|5|i. 
This limit is equivalent to keeping only the terms of order N in the 1/A expansion. We observe from 
BCS expressions (j5J) that the expansion in 1/A converges for A > 1/ir. In the opposite case of one 
pair and two levels, 2M = N = 2, the ground state energy can be computed exactly by e.g. solving 
Eqs. (J3J) with the result 

E 2 gT = -d(\ + VTTV) (22) 

In this case the expansion of the ground state energy ((22]) in 1/A converges for A > 1. In general, we 
believe that strong coupling expansion yields convergent rather than asymptotic series with the radius 
of convergence between A c i ~ 1 and A C 2 ~ 1/ir. 

Bellow in this Section we develop an efficient algorithm for calculating the low energy spectrum 
to any order in 1/A. While the pseudospin representation detailed in the previous Section provides a 
simple and intuitive description of the strong coupling limit, the usual perturbation theory becomes 
unmanageable beyond the first two orders in 1/A. An approach based on Bethe's Ansatz equations, 
on the other hand, turns out to be well suited for the purposes of systematic expansion. 
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3.1 The ground state 



Here we expand the ground state energy in 1/A. Richardson's equations Q lead to recurrence relations 
for the coefficients of the expansion. From these relations the ground state energy can be computed 
to any reasonably high order in 1/A, e.g. we write down the energy up to 1/A 7 . As it was mentioned 
above we take the number of electrons to be even and consider only the case when Xd > J. As we have 
seen in the previous Section, this inequality ensures that in the ground state all levels are unblocked 
and all electrons are paired, i.e Richardson's equations (jSJ) should be solved at 

m = M n = N 



We begin by introducing a convenient set of variables 

N M i 

s p = Y.^ k ) p o^Epp (23) 

k=l i=l 

Variables a p can be expanded into series in the inverse coupling constant A. 

oo 

°v = 12 a >~ k ~ P ( 24 ) 

k=0 

Next, we rewrite Richardson's Eqs. Q in a form more suitable for our purpose. We divide the equation 
for Ei by Ef with p > — 1 and add all M equations for each p. Expanding 1/(1 — 2e k jE^) in le^jEi 
and using an identity 

2 f 1 1 A p 

i>3 1 3 X^i / k=l 

we obtain 



M 

E gT (M,N,s p )=J2Ei 



i=l 



oo / j+1 

J-k+1 



(25) 



—M(N -M+ l)Xd - s k<?k = —M(N -M+ l)Xd - dJ2 [12 s k4 ] X~ j 

k=l j=0 \k=l 



a v 

~ w - 12 °v-k+iOk = (N- p)a p+ i + 22 s i a 3+P+i V > (26) 

k=l j=l 

Now plugging a p = J2 < kLo a p^~ k ~ p the last equation and setting the coefficient at A~ h ~ p_1 to 

h p h h 
P 

d 



zero, we obtain 



7 + EE + £ s k a h p - k k+1 — —(N — p)a h p+1 (27) 

k=l s=0 k=l 



Note that from <7o = M it follows a° = M and = for k > 1. The values of a,Q serve as boundary 
conditions for recurrence relations f)27j) . Note also that according to (|27j) the coefficients a p h do not 
depend on A as expected from their definition (J2U). Coefficients a p determine a p for the ground state 
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to the lowest nonvanishing order in 1/ A and therefore can be expressed in terms of zeroes of Laguerre 
polynomial (|21|) with r = M. Using (|21|). we obtain 



al<P = {-iy—\ogL- A t N {x) 



x=0 



According to Eq. (|23j) in order to determine the ground state energy to order 1/A J one has to calculate 
the first j — p + 2 coefficients in the expansion of a p . To do this, we first compute a° for p < j + 1, 
then dp for p < j, then a 2 for p < j — 1, etc. In other words, we start from a? element of matrix ap 1 and 
use recurrence relations (|2"Tj) to move down the first column of this matrix until a® +1 , then to move 
down the second column from a\ to a] etc. 

While we were not able to express a p in terms of p and k explicitly, the above procedure allows for 
an efficient calculation, e.g. using Mathematica, of the ground state energy to any given order. For 
example, the ground state energy to order 1/A 2 is 

E gI (M, N, „) = -M(N - M + 1) « + £# - U - d) M{N ~ M » 



N V N J N 2 (N - l)Xd 

sl\ M(N — M)(N — 2M) f Sl s 2 \ M(N — M)(N — 2M) 

NJ Sl N 4 {N - l){Xd) 2 + V 3 ~ ~WJ N 3 (N — 1)(N — 2)(Xd) 2 

(28) 

From Eq. ()28|1 one can make several observations. 

1. For N = M the first two terms give the exact energy. This is seen by noting that N = M 
means that all levels are doubly occupied, i.e. there is only one state. Averaging Hamiltonian (0) over 
this state gives the exact energy of the system which turns out to be equal to the first two terms in 
(}2"%j) . Therefore, the remaining terms in the 1/A series for the ground state energy are proportional to 
N — M. 

2. When iV = 2M, all terms with even nonzero powers of 1/A vanish. This can be demonstrated, 
e.g., by writing the kinetic term in the BCS Hamiltonian (J2J) in terms of pseudospin operators © 

N 

H(X = 0) = 5>^f + = H + -1 (29) 
i=i / A 

and noting that iV = 2M correspond to zero ^-projection of the total pseudospin. In this case, by 
Wigner-Eckart's theorem K? has nonzero matrix elements only for transitions K — > iT± 1, while 
matrix elements for transitions K — > K are equal to zero. The terms with even nonzero powers of 
1/A vanish because they contain at least one matrix element of Hq from (|2T)|) between states with 
the same K. These terms are therefore proportional to iV — 2M . Even terms also vanish when are 
distributed symmetrically with respect to zero. Hence, they reflect an asymmetry in the distribution of 
€i. For example, the ground state energy for iV = 2M and equidistant single electron levels distributed 
symmetrically between ±D = ±(m — l/2)d is 



E 2m = -Dm 



^2m + 2 2m + 1 16m 2 + 22m + 7 
2m- 1 + 3(2m- 1)A ~ 180(2m - 1) 2 A 3 + 

(30) 

128m 3 + 380m 2 + 344m + 93 _ /x7N 

- + 0(1/A 7 ) 



7560(2m - 1) 3 A 5 



9 



One can check that in the limit m — > oo this expression reproduces the BCS result (j5J) for the ground 
state energy up to terms of order 1/A 7 , while for m = 1 we recover (|22j). Note also that the case p = N 
in Eq. ()27|) does not seem to be problematic as at p = N the factors of 1/ (N — p) in Eqs. (j23 IHUjl are 
always compensated by a factor of (AT — p) in the numerator of the corresponding term. 

3. Richardson's equations Q remain invariant if single electron levels e k are shifted by 5 and 
parameters Ei are shifted by 25. The total energy E = Y$±Li Ei then shifts by 2Mb. Note that this 
shift is entirely contained in the second term of expansion (J28)) . Thus, the remaining combinations of 
s k at each power of 1/ A are "shiftless" . For example, 

2 2 

s 2 - ^ - s 2 + 26 Sl + N6 2 - (si + 2N6 Sl + N 2 5 2 )/N = s 2 - ^ 
3.2 Excited states 

Let us now expand energies of low-lying excitations in 1/A. These expansions turn out to be analogous 
to that for the ground state energy. We begin with the excitations that conserve the number of pairs 
and then turn to the simpler case of pair-breaking excitations. 

It was demonstrated in Section 2 that for Xd > J lowest pair-conserving excitations correspond to 
total pseudospin K = N/2 — 1 and total spin S = 0, where iV is the total number of single particle 
levels. The number of such states according to degeneracy formula (fT2"j) is iV — 1 and their energy is 
—XdK(K + 1) according to (jOJ). We also know from Section 2 that for these states one of parameters 
Ei (say Em) remains finite as A — ► oo, while all others diverge in this limit. 

To distinguish Em from the rest of parameters Ei, we denote it by rj. Richardson's equations (J3J) 
read 

i M-l 9 N 1 9 

f = J2 — KM (31) 

Xd . =1 Ei — Ej k=1 Ei — 2efc Ei — rj 

1 M-l 9 Af i 

- -L - V = V — !— i = M (32) 

Expanding the LHS of Eqs. (|31|) in 2e k / Ei and rj/Ei and performing the same manipulations that 
lead to Eqs. (J251 I2SJ) for the ground state, we obtain 

M-l oo 

Y,Ei= ~( M - 1)(JV - M)Ad - - 2t/>* (33) 

i=i fe=i 

a, p 00 
- - X] Vp-k+iVk = {N -p- 2)a p+1 + X( s i - 2v J )° r j+ P +i V > (34) 
A " ik=i j=i 

where now a p = ^fii 1 V^f- We see that replacements 

M -»• M - 1 iV -> iV - 2 s p -> s p - 2rf (35) 

transform Eqs. ()331 13*1^) into Eqs. (|23l l2T?|) for the ground state. Thus, energies of first N — 1 excited 
states are 

M-l 

^pair =J2 E i + E M = E gT (M - 1, N - 2, s p - 2rf) + V (36) 
/=i 
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Let us also rewrite Eq. (|32|) for rj as 



N x j 

E = " w " 2ff i " 2 W2 - 2 V 2 a 3 - . . . (37) 

One can see (by e.g. sketching the LHS of Eq. (|37|) ) that this equation has iV — 1 roots with the kth 
root lying between 2e k and 2efc+i . To the lowest order in 1/A this equation reads 

N 1 

Y = (38) 

Eqs. (JBH) and (|3Tjl are to be solved iteratively order by order in 1/A. The procedure is similar to 
that for the ground state, e.g., recurrent relations analogous to (J2"7|) can also be derived. The only 
difference is that the coefficients at powers of 1/A now depend also on i] Q , which has to be obtained 
from (|38J). For example, the excitation energies (JBljj) to the first two orders in 1/A are 

^pair = -(M -l)(N- M)Xd + ^ " " ^ + Vo (39) 

^pair - % = NXd + 770(1 - 2f) (40) 

where 

f = (M- 1)/(JV - 1) w M/AT (41) 

is the filling ratio. 

Energies of higher excitations can be computed in the same way by solving 2, 3, 4, ... coupled 
equations of the type of (|3~%j). For instance, energies of the next group of excited levels to the first two 
orders in 1/A are determined by solutions of the system 

t 1 



k =i Vi ~ 2e fc fjt - 7]2 



E 



\ r} 2 - 2e k r)! - r] 2 

Now let us consider pair-breaking excitations. For Xd > J low energy excitations of this sort 
correspond to breaking a single pair of electrons thereby decreasing the number of pairs by 1 and the 
number of unblocked levels by 2. Let the single electron levels occupied by two unpaired electrons have 
energies e a and e^. Since the lowest energy is achieved by having the unpaired electrons in a triplet 
state (recall that J > corresponds to the ferromagnetic exchange), the energy of lowest pair-breaking 
excitations according to (jH) is 

^spin = e a + e b -2J + E gT (N - 2, M - 1, s p - (2e a ) p - (2e b ) p ) (42) 

Note that, unlike i] in (pjj) . single electron energies e a and e b do not depend on A. Therefore, to 
compute the energy of pair-breaking excitations we need only recursion relations ()27|) for the ground 



11 



state with N' = N — 2, M' = M — 1, and s' = s p — e p a — . In particular, to the first two orders in 
1/A we get from (J5BJ) 

p (m i\(*r A/f\ , (si-2e a -2e b )(M- 1) 

E spin = -( M - ^(^ - M )# + + e a + e b - 2J 

^spin - E gr « ^Ad + (e« + - 2/) - 2 J (43) 

It is instructive to compare the above results with the BCS theory [§]. For this purpose let us write 
down the energies of the pair-conserving excitations for large M and iV up to the order 1/A. 

^pair - ^gr « + Vk (l - 2/) + rf j^'^ (44) 

where -D = Nd, f is the filling ratio 1)41 J), and % is the fcth root of Eq. ()37)1 . In deriving the above 
equation from (j4*U|) and (j2~H|) we shifted the single electron levels so that q = (j2iLi e ij /N = 0. In 
BCS theory (i.e. in the limit N,M —>■ 00) pair-conserving excitation energies are |XU] 



2 V ( e k - v) 2 + A 2 (45) 

where // is the chemical potential and A is the gap. In the strong coupling regime both \i and A are 
of order A. Expanding the square root in expression (J4*5j) in small up to e|, we see that (|43|) and 
coincide to this order if we identify 



Vk = 2e k A = DXy4f(l - f) fi = (2f-l)D\ 

The first of these equations follows from (p7f|) in the limit of large N, while the remaining two can 
be derived from the BCS equation for the gap and chemical potential (see e.g. JHI)- Similarly one 
can check that pair-breaking excitations ()4*2^) correspond to two Bogoliubov quasi-particles with total 
energy 



e a - f ,f + A^ + ^(e b - f xf + A^ (46) 

Note that in the BCS limit the difference between pair-breaking and pair-conserving excitations 
disappears and expression ()43j) simply corresponds to two quasi-particles in a singlet state each having 



the energy yj (e k - /i) 2 + A 2 . 

We have seen in Section 2 (see also Fig. 1) that in the strong coupling limit many-particle energy 
levels of the BCS Hamiltonian (J2J) coalesce into narrow well separated bands. Expression ()4*4*|) can be 
used to estimate the ratio of the width of the first band, Wi, to the single particle bandwidth D = Nd. 

Wi f(l - f) 

^ « 2(1 - 2f) + (47) 

where W\ is the width. Note that at half filling, / = 1/2, the width of the first band goes to zero as 
A — > 00. In general, it follows from Wigner-Eckart's theorem (23] (see the discussion in item 2 under 
the ground state formula (J2EJ)) that at half filling widths of higher bands also vanish as A — > 00. 

According to the BCS equations for the excitation energies and (j4*Uj) the gaps A g pj n = 



E„„:„ — Ep 



and A 



spin S r J m j„ ^pair " Impair -^g r 



for the two types of excitations coincide in the 
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thermodynamical limit. We have also seen in Section 2 (see the discussion bellow degeneracy formula 
([12))) that when < J < Xd and J/ (Xd) remains finite as A — > oo, spin-1 excitations have lower energy 
as compared to pair-conserving excitations . If, however, J ~ d or smaller, keeping J to the lowest 
order in 1/A in excitation energy ()42j) is not justified. In this case the two gaps are the same to this 
order. Therefore, it is interesting to set J = and evaluate the gaps to the next nonzero order. 
Depending on the filling ratio / (see Eq. ([41))) we can distinguish two different cases. 

1. / ^ 1/2. Lowest lying excitations correspond to smallest or largest possible values of i] and 
e a + 6;, depending on the sign of (1 — 2f). To determine the maximal and minimal r] Q , note that the kth 
root of Eq. (J3~5|) lies between 2e^ and 2e&+i. If iV is large and — e^+i — > as iV — > oo, the smallest 
and largest solutions of ()38j) are rj^ im pa 2e\ and n™ ax pa 2e n respectively. We have from (j4"Ul l4"3*)) 

^spin" A pair = ci|l-2/| >0 (48) 

where we have used e n — e n _i ~ £2 — ei ~ d and d is the mean level spacing. 

2. / = 1/2. To the first two orders in 1/A: A g pj n — Ap a j r = 0. In the next order we obtain from 

(EH1ESI andi2) 

A • - A • = ^ f!« f£6 (49) 

^pair ^spm 2 jVAd 1 J 

where we shifted single electron levels so that e$ = fe^li e ij /N = 0. We show in the Appendix using 
Eq. ()38|) for r] that the minimal value of t]q is always smaller than that of 2(e^ + e\). Therefore, 
^spin > ^pair- 

Thus, at J = the pair-breaking excitations always have a larger gap in the strong coupling limit. 
Note that for A = the situation is opposite as it always costs less energy to move one of the two 
electrons on the highest occupied single electron levels to the next available level. Since according to 
BCS expression (J5J the energy gap in the strong coupling limit is 2A pa 2DX = NXd, we see from (149)) 
that at the half-filling Ap a j r — A g pj n pa d 2 /A, i.e. the difference between the two gaps vanishes in 
the thermodynamical limit. 



4 Conclusion 

We determined the spectrum of the Universal Hamiltonian {1} in the strong superconducting coupling 
(A > 1) limit (0 E21 EI} and developed a systematic expansion in 1/A around this limit (|2*7l I2~%1 
I3TH |4"2*|) for the ground state and low-lying excitation energies. We detailed an algorithm by which 
these energies can be explicitly evaluated up to arbitrary high order in 1/A and estimated that the 
expansion converges for A > A c where A c lies between A c i ~ 1 and A C 2 ~ 1/ir. Technically, this 
expansion is based on the existence of the exact solution jB] of the BCS Hamiltonian (J2J). We found 
that in the strong coupling limit Richardson parameters are zeroes of appropriate Laguerre polynomials 
(121)) and analyzed their behavior at large enough but finite A . 

We found that it is important to distinguish between two types of excitations in the problem: those 
that conserve the total number of paired electrons and those that do not. We determined the energy 
gaps for both types and found that at zero spin-exchange constant, J = 0, in contrast to the weak 
superconducting coupling limit, the gap for pair-breaking excitations is always larger ([48 ) I49| ) . 

We believe there are two physically motivated questions within the scope of validity (see [3]) of 
the Universal Hamiltonian (0) that still need further clarification. The first problem is to develop 



13 



a quantitative description of the crossover between a perturbed Fermi gas and the region of strong 
superconducting correlations (see [T3] and the discussion in the beginning of Sections 2 and 3). The 
second problem is to study analytically the interplay between superconducting correlations and spin- 
exchange (see e.g. 



5 Appendix 



We show here using Eq. ()38|) for 770 that the minimal value of rj 2 is always smaller than that of 2(e 2 + elj, 
i.e. 



xl < 2(a 2 + b 2 ) (50) 

where Xq is the smallest in absolute value solution of (|38j) . a and b are the two smallest in absolute 
value single electron levels 6j, and \a\ < \b\. Indeed, consider a function 

N 1 

»W-£jraS (51) 

To prove f)5U|) we need to show that g(x) has a zero on the interval (— c, c), where 



c = ^2(a 2 + 6 2 ) 

For N = 2 there is only one zero, Xq = E\ + £2, and (|5U|) clearly holds. Consider N > 2. 

8(x)k 




Figure 2: A schematic plot of the function g(x) = J2k=i x -2e k on ^ ne m t erva l from — c to c, where 
c = \j2(a 2 + 6 2 ), a and 6 are the two smallest in absolute value single electron levels e^, and \a\ < \b\. 
Note that since 2\b\ > c there is only one pole on this interval. In the vicinity of 2a we have g(x) ~ 
l/(x — 2a) and therefore g(x) is positive on the immediate right of x = 2a and negative on the left. 

First, note that g(x) has a single pole at x = a on this interval from — c to c, and g{a+) > 0, while 
g(a—) < 0. Hence, there is a zero between c and — c iff either g(c) < or g(— c) > 0. To show that 
this is the case it is sufficient to demonstrate that g(c) — g(— c) < 0. We have 

N 2c 2c 

g( C )-g(-c) = Y,^r^2= E ^3^2 
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which is indeed negative since c 2 < 4e 2 for all £j except = a, b. 
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